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If the cosmological QCD phase transition is strongly first order and lasts sufficiently long, it gener- 
ates a background of gravitational waves which may be detected via pulsar timing experiments. We 
estimate the amplitude and the spectral shape of such a background and we discuss its detectability 
prospects. 
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I. INTRODUCTION 

Gravitational waves (GWs) are space-time fluctuations 
that propagate at the speed of light through empty space; 
they were predicted by Einstein in 1916 (y. Because of 
the weakness of the gravitational interaction, GWs could 
provide information about astrophysics and cosmology 
from regions and epochs of the Universe from which elec- 
tromagnetic radiation cannot propagate freely. For the 
same reason, however, GWs have thus far eluded direct 
detection, despite considerable efforts. 

Advanced configurations of existing ground-based in- 
terferometers such as LIGO [2| and VIRGO [3[ are ex- 
pected to detect GWs in the next years. Terrestrial in- 
terferometers have best sensitivity at a frequency / ~ 
100 Hz, and are severely limited by seismic noise below 
a few Hertz. GWs with significantly lower frequency, 
/ ~ 10~ 9 Hz, arc also expected to be detected by pulsar 
timing experiments in the next decade [J, Q ■ A world- 
wide collaboration of astronomers, the International Pul- 
sar Timing Array (IPTA) project jj], has been formed 
with the goal of detecting nano-Hz GWs using millisec- 
ond pulsars. Millisecond pulsars are rapidly rotating, 
highly magnetized neutron stars which emit a beam of 
electromagnetic radiation that sweeps over the earth once 
per rotation. They constitute extremely accurate clocks 
that could be used to detect GWs. Candidates for the 
generation of a GW background in the nano-Hz band 
are super-massive black hole binary mergers |6l— Till] and 
cosmic strings [I2I - FI6J . 

In this paper we study another potential candidate: 
the cosmological QCD phase transition, which is believed 
to have taken place when the Universe had a temperature 
of T* ~ 100 MeV. A GW background can be generated 
by the QCD phase transition if it is first order, and the 
characteristic frequency of this background falls in the 
frequency band of pulsar timing experiments. We show 
here that, if the phase transition is sufficiently strong and 
lasts for a sufficiently long time, the GWs produced can 
be observed in future pulsar timing experiments. This 
possibility has been discussed for the first time by Wit- 



ten in Ref. [17|. Here we present accurate predictions 
for the spectrum of the emitted gravitational radiation 
as a function of the phase transition parameters like its 
temperature, strength and duration. 

In the context of standard cosmology and QCD, the 
cosmological QCD phase transition is not even second 
order but a cross-over, and we do not expect it to gener- 
ate GWs. However, if the neutrino chemical potential is 
sufficiently large (still well within the bounds allowe d by 
big bang nucleosynthesis) , it can become first order [18J] . 
Furthermore, if a sterile neutrino is the dark matter, we 
do expect a large neutrino chemical potential [19| . 

Thus, pulsar timing experiments could open a new cos- 
mological window: the detection of a stochastic back- 
ground of GWs could help to determine whether the 
QCD phase transition is first order. The amplitude and 
peak frequency of the spectrum are also sensitive to the 
expansion rate of the Universe during this phase transi- 
tion [20j, which is currently unconstrained. 

In the next section we provide estimates for the GW 
spectrum by a first order QCD phase transition. In Sec- 
tion IIIII we compare our results with current and ex- 
pected sensitivities of pulsar timing arrays. We conclude 
in Section IIVI Throughout we use the metric signature 
(—,+,+,+) and conformal time so that the Friedmann 
metric is given by ds 2 = a 2 (t) (— dt 2 + dx 2 ). The confor- 
mal Hubble parameter is denoted by rl = a/ a = Ha. An 
over-dot denotes the derivative w.r.t conformal time t. 



II. GRAVITATIONAL WAVES FROM A FIRST 
ORDER QCD PHASE TRANSITION 

Very violent processes in the early universe can lead to 
the generation of GWs. One example of such violent pro- 
cesses are first order phase transitions [13, I2ll - l24j , which 
can lead to GW production via the collision of bubbles of 
the true vacuum [25l - l33J | and via the turbulence and mag- 
netic fields they can induce in the cosmic plasma [34144 1| . 

The GWs generated by a source are determined by the 
linearized Einstein equation for tensor perturbations in a 



Fricdmann background [42 



h± + 2-h± + k 2 h± 



87rGa 2 pn d 



(1) 



Here p is the background cosmological energy density, a is 
the scale factor, H± are the two tensor helicity modes of 
the (dimensionless) anisotropic stress which is the source 
of GWs, and h± are the helicity modes of the GW. If 
an anisotropic stress is generated at some time £* in the 
radiation dominated era, from a source with relative en- 
ergy density Qg* — Ps*l P*, wc expect II to be at best 
of the order of fig*. The GW energy density (from two 
polarizations ± which contribute equally) is given by 



PGw(t) = 



(h + (x,t)h + (x,t)) 
8nGa 2 (t) 



(2) 



Due to statistical homogeneity, pgw is independent of 
the position. The GW energy spectrum per logarithmic 
unit of frequency ^^ \ k \ is defined by 



Pgw 



dk dpGw(k,t) 
k dlog(fc) 



Detailed semi-analytical and numerical calculations have 
been performed in the past in order to calculate the GW 
energy spectrum from first order phase transitions |25l - 
1-411 ] . In this paper we simply use analytic fits to the most 
recent results, as presented in the following. 

Concerning the GW signal from bubble collisions, we 
use the shape of the spectrum proposed in Ref. [33| s but 
rescale the amplitude to agree with the numerical result 
of Ref. [3jJ. As a result, the GW energy density emitted 
by this source is well approximated by 



Bubble collisions: 
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Here f2 r o denotes the radiation energy density today, h = 
i?o/(100km/s/Mpc) is the present Hubble parameter in 
units of lOOkm/s/Mpc, /3 _1 is the duration of the phase 
transition, v is the expansion velocity of the bubbles, 
and k is the co-moving wave-number or frequency of the 
GW. The GW spectrum is proportional to the relative 
energy density in the source f2|„, to the ratio between 
the duration of the phase transition and the Hubble time 
(H*//3) 2 , and to the bubble velocity v 3 HJ. 

The QCD phase transition is expected to happen 
at the temperature T* ~ 100 McV, when the kinetic 
and magnetic Reynolds numbers of the cosmic fluid 
are very large [4l| . The bubbles which rapidly expand 
and collide are therefore expected to generate magneto- 
hydrodynamical (MHD) turbulence in the cosmic fluid. 
The kinetic energy of the turbulent motions and the mag- 
netic fields sustained by the MHD turbulence also induce 



GWs: Ref. [4l| presents the latest semi-analytical calcu- 
lation of the GW spectrum from MHD turbulence. There 
are two important differences with respect to the GW sig- 
nal from bubbles. First, turbulence lasts beyond the du- 
ration of the phase transition: this leads to an enhance- 



ment of the signal on large (super-horizon) scales [41 



Second, the time correlation properties of the anisotropic 
stress source are different. For bubble collisions, the 
source is totally coherent (see [33|, [41( ) , while for MHD 
turbulence the source is coherent only over about one 
characteristic wavelength [41|. This leads to a difference 
in the peak position of the GWs from the two sources: 
while the signal from bubble collisions peaks at k p ~ j3 t 
the inverse duration of the phase transition, the peak of 
the MHD signal is related to the bubble size: the peak 
wavelength becomes therefore A p ~ i?* ~ v/f3. The anal- 
ysis of Ref. [4jJ finds a peak at about k p ~ ir 2 /3/v. We 
can fit the GW spectrum obtained in 4l| by the following 
formula: 



MHD turbulence: 
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For large scales, k <C k p , both spectra in Eqs. ^ 
and ([5]) increase as k 3 : this behavior is simply due to 
causality [33, [43| . Since the anisotropic stresses are gen- 
erated by a causal process, their spectrum is white noise 
at scales larger than the typical correlation scale of the 
source, which corresponds to the bubble size. The white 
noise spectrum is inherited by the GWs: (\h\ 2 ) oc const. , 
so that the GW energy density scales simply with the 
phase space volume k 3 . The behavior on small scales, 
k ^$> k p , depends on the source power spectrum and on 
the unequal time correlation properties of the source, see 
Ref. (4J|. In particular, the result of Eq. ^ resides on 
the assumption that the bubbles are infinitely thin: this 
assumption holds if the bubbles propagate as detonations 
and causes the fc -1 slope at high wave-numbers [3l[[33|. 
On the other hand, the fc -5 / 3 decay of Eq. (01 is a con- 
sequence of the Kolmogorov type spectrum assumed for 
the MHD turbulent motions at high wave-numbers. In 
addition, the slope of the MHD signal changes at sub- 
horizon scales, H* < k < k p , from fc 3 to k 2 due to the 
long duration of the source (c.f. [4l]|). 

From the the above formulas for the GW spectra we 
see that the basic ingredients which determine the peak 
position and amplitude are simply the fractional energy 
density of the source fls*, the duration of the phase 
transition /3 _1 , and the bubble velocity v (besides ob- 
viously the temperature at which the phase transition 
occurs T*, which is parameterized by the Hubble scale 
%*). Q.s* and v are related, in a way which depends on 
the characteristics of the phase transition: for example 
its strength, the properties of the bubble expansion, the 



interactions of the fluid particles with the field which is 
undergoing the transition, and so on. In early works on 
GWs from bubble collisions, it has been assumed that 
the bubbles expand as a Jouguet detonation, because in 
this case the above parameters can be calculated quite 
straightforwardly [27|. However, the recent analysis of 
}4a] demonstrates that there is no particularly justified 
reason for this assumption, and other kinds of solutions 
for the bubble expansion are possible, such as deflagra- 
tions, runway solutions and hybrids. Ref. [45| presents 
a model-independent description of the different regimes 
characterizing the bubble expansion, including the effect 
of friction due to the interaction of the bubble wall with 
the fluid particles. From a given particle physics model, 
one can in principle evaluate the friction parameter r\ and 
the strength of the phase transition a = p vac / Prad- Once 
these two quantities are known, Ref. [45[ provides a way 
to determine the bubble wall velocity v and the fraction 
of vacuum energy density which goes into kinetic energy 
of the bubble walls, K = Pkin/Pvac- In terms of the param- 
eter k, the fractional source energy density for bubbles 
becomes 



O 



(B) 
S* 



a + 1 



In the case of MHD turbulence, one further has to convert 
a part of the bubble wall kinetic energy to turbulence and 
magnetic fields. The efficiency of this conversion is not 
straightforward to estimate. Ref. [45| provides a relation 
between the bubble wall velocity and the fluid velocity 
at the bubble wall position vf. in the most optimistic 
case, one can argue that the overall kinetic energy in 
turbulence is simply determined by this fluid velocity. In 



this case one would have Q~ ' 



vV2. 



In the absence of a way to determine a and rj from 
a given particle physics model, in the present analy- 
sis we have decided to keep the parameters completely 
model- independent. We make the simple assumption of 
equipartition, namely we assume the same energy den- 
sity in colliding bubble walls and in MHD turbulence, 
n s J = Q s ^ . This is more a reflection of our ig- 
norance of how this energy density will be distributed 
than a well justified assumption; nevertheless, it seems 
to be a reasonable expectation and in the absence of a 
model for the phase transition it is the most straightfor- 
ward assumption. We also assume a strongly first order 
phase transition, which induces super-sonic bubble ve- 
locities, v > c s . We set the temperature of the QCD 
phase transition at T* = 100 MeV. The other parame- 
ter relevant for the GW spectra is the duration of the 
phase transition, parametrised by j3. In the electroweak 
case, this is usually taken to be (1 - 10)% of a Hubble 
time: /? = (10 — 100)%*. This value is based on the 
estimate given in Ref. [2l|, which shows that /? is re- 
lated to the temperature of the phase transition through 
fl/H* ~ 41n(mpi/T»), for a phase transition nucleated 
via thermal fluctuations. In the absence of a precise 
model for the QCD phase transition, in this section we 



have decided to take /3 = 10 %*, which is more favorable 
for observations with pulsar timing experiments. The 
analysis of |2l[ demonstrates that models of phase tran- 
sitions with small values of /3/'H* may be rather excep- 
tional, but cannot be ruled out by general arguments. 
The important point is that j3/'H* must be larger than 
unity, otherwise the phase transition is not fast with re- 
spect to the universe expansion and our assumptions no 
longer hold. 




FIG. 1: The GW spectra from bubble collisions (black, solid) 
and from MHD turbulence (red, dashed) are shown for dif- 
ferent values of f2s* =0.1 and v = 0.7 (top panel) and 
Qs* = 0.03 and v — 0.57 ~ c s (bottom panel). We set 
P = 1QH, and T, = 100 McV throughout. 

In Fig. Q] we show the GW spectrum for both, bub- 
bles and MHD turbulence for two different choices of 
the parameters Vis* an d v. The MHD turbulence sig- 
nal dominates almost in the entire frequency range. At 
large scales, it is slightly higher due to the long duration 
of the turbulent source with respect to bubble collisions 
[4l| . As already mentioned, the long duration of the 
source also causes the slope of the MHD signal to change 
at sub- horizon scales from k to k 2 : consequently, for 
(3 > k > H*, i.e. 0.1 < k/f3 < 1, the bubble colli- 
sion signal prevails. This is valid up to the peak of the 
bubble collision signal, which arises before the turbulent 
one: at k/j3 ~ 1, corresponding to the inverse charac- 
teristic time of the source, while the turbulent spectrum 
peaks at k/j3 ~ 7r 2 /t>, corresponding to the inverse char- 
acteristic scale of the source. This causes the turbulent 
signal to dominate at interesting frequencies, since the 
total spectrum continues to raise after k//3 ~ 1 (only if 
the energy in turbulence is about one order of magnitude 
smaller than the one in bubble collisions, the collision 



signal will dominate: however, this seems somewhat un- 
natural given the extremely high Reynolds number of the 
primordial fluid, and we discard this possibility in this 
work) H. 
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FIG. 2: The GW signal from bubble collisions and MHD tur- 
bulence for Q. s , = 0.1 and v = 0.7. We choose /3 = WH*. 
The signal is dominated by the contribution from MHD tur- 
bulence. The bubble collision peak causes the hump on the 
left of the true peak of the spectrum. 

In Fig. [5] we show the total signal for the more opti- 
mistic case, Qs* = 0.1 and v = 0.7. The peak frequency 
of the total GW spectrum corresponds to the MHD tur- 
bulence peak: k/(3 ~ 7r 2 /u, and depends on the choice 
/3 = lOK*. From / = k/(2%) one obtains gHH] 



f P s 1.7 ■ 10- 



J3_ 
H* 



m 



r. 



100 MeV 



Hz 



(5) 



where g* is the number of effective relativistic degrees of 
freedom at the temperature T* . Withv = 0.7, (3 = 10 H*, 
<7* = 10 and T* = 100 MeV the peak frequency becomes 
f p ~ 2.5 • 10~ 7 Hz. 



III. THE PULSAR TIMING ARRAY 

Neutron stars can emit powerful beams of electromag- 
netic waves from their magnetic poles. As the stars ro- 
tate the beams sweep through space like the beacon of a 
lighthouse. If the Earth lies within the sweep of a neu- 
tron star's beams, the star is observed as a point source 
in space emitting short, rapid pulses of electromagnetic 
waves, and is referred to as a pulsar. 

The electromagnetic pulses we observe arrive at a very 
steady rate due to the enormous moment of inertia of 
neutron stars. The idea to use these stable clocks to 
detect GWs was first put forward in the late 1970s [471 - 
|49I |. Fluctuations in the time of arrival of pulses, after all 
known effects are subtracted, could be due to the pres- 
ence of GWs. Recently pulsar timing precision has im- 
proved dramatically. Jenet and collaborators [50J have 
shown that the presence of nano-Hertz GWs could be 
detected using a pulsar timing array (PTA) consisting 



of 20 pulsars with timing precisions of 100 nanoseconds 
over a period of 5 to 10 years (see also [J, Q for more re- 
cent PTA sensitivity estimates) . Pulsar timing arrays are 
most sensitive in the band 10 ~ 9 Hz < / < 10~ 7 Hz. The 
lower limit in frequency is given by the duration of the ex- 
periment (~ 10 yr.) and the upper limit by the sampling 
theorem, i.e. the time between observations (~ 1 month). 
The spike in the sensitivity at / = 0.3 x 10 _7 Hz seen in 
Fig. [3] is the frequency of the earth's rotation around the 
sun which cannot be disentangled from a GW with the 
same frequency. 



Current NANOGrav sensitivity 

PTA 2020 




FIG. 3: Comparison of the GW spectrum h 2 Q(f) with cur- 
rent NANOGrav pulsar timing array sensitivity and expected 
sensitivity of pulsar timing experiments in 2020 [5| . We have 
used h = 0.73, Q r0 = 8.5 x 10~ 5 , tt s * = 0.1 and v = 0.7. We 
plot the GW spectra for the values H*//3 = 1, 0.5, 0.2, 0.1 
(dashed lines from top to bottom). For 'H,/j3 ~ 1, the back- 
ground of GWs can just be detected in present pulsar timing 
experiments, while for 0.1 < %*//3 it can be detected by the 
planned array IPTA2020 (very high values of W*//3 ~ 1 are 
difficult to accommodate in the case of a thermally nucle- 
ated phase transition, c.f. discussion in the text). We also 
show the LISA sensitivity [52l . [53J . Unfortunately, LISA will 
not be able to detect a signal from a first order QCD phase 
transition (the EW phase transition is more promising in this 
respect [25M4ll . lil. |4^| ) . 

The North American Nanohertz Observatory for Grav- 
itational Waves (NANOGrav) j51|, a collaboration of as- 
tronomers, has created a pulsar timing array-a galactic 
scale GW observatory using about 20 pulsars. It is a 
section of the IPTA, an international collaboration in- 
volving similar organizations of European and Australian 
astronomers. The current NANOGrav pulsar timing ar- 
ray sensitivity is shown in Fig. El together with the GW 
spectra we expect from the QCD phase transition as a 
function of frequency 



h 2 n GW {f) = k 



dVL 



GW 



d log k ' 



(6) 



for H*/l3 = 1 (top dashed line), H*/f3 = 0.5 (upper- 
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FIG. 4: Same as Fig.|3]but comparing the characteristic strain 
h c as given by Eq. (0). 



middle dashed line), H*/j3 = 0.2 (lower-middle dashed 
line), H*/f3 = 0.1 (bottom dashed line). We have used 
h = 0.73, f2 r o = 8.5 x 10~ 5 (which includes photons and 
neutrinos), £1$* = 0.1 and v = 0.7. We have taken the 
frequency to be (c.f. Eq. (J5J) 



/ = 1.7x 10" a £- 



l n 



i.e. chosen T* = 100 MeV, g* = 10, and k/0 varying 
between 10~ 4 and 10 4 (like in Figs Q] and 0). The signal is 
compared with the current sensitivity of the NANOGrav 
pulsar timing array, and the expected sensitivity of the 
IPTA pulsar timing array in 2020 j5|. For values of 0.1 < 
W-*/P ^ lj the background of GWs would be detected 
with future pulsar timing array sensitivities. The value of 
Tt*/ fi must certainly be smaller than unity for the phase 
transition to be fast with respect to the Hubble time and 
for our approximations to apply. In most cases, if the 
phase transition happens at a temperature much smaller 
than the Planck temperature, H*//3 is of the order of 
0.01; however, higher values of this parameter cannot be 
excluded, and we adopt them here since they are more 
promising for detection [21J . 

Fig. |3] also shows the sensitivity of the planned Laser 
Interferometer Space Antenna (LISA) [52| assuming that 
some of the confusion noise from white dwarf binaries 
can be subtracted out [53( . LISA will not be able to de- 
tect the GW signature of a first order QCD phase tran- 
sition: in order to be detectable by LISA, the GW spec- 
trum must peak at higher frequency and consequently the 
phase transition must occur at higher temperature. LISA 
can in principle detect GWs from a strongly first order 



EW phase transition at T* ~ 100 GeV [25H411 |4J, IM |5J| . 
A related quantity often used in the pulsar timing com- 
munity is the (dimensionless) characteristic strain, de- 
fined by [H 



hl(f) 



3H§ 

2ir 2 
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(7) 



In Fig. 2] we show the same data as in Fig. [3] but in terms 
of the characteristic strain h c . 



IV. CONCLUSION AND OUTLOOK 

A stochastic background of GWs from the QCD phase 
transition could be detected by pulsar timing experi- 
ments if the transition is strongly first order and lasts 
sufficiently long with respect to the Hubble time. In 
standard cosmology the QCD phase transition is not even 
second order, but simply a cross-over and in this case we 
do not expect it to generate GWs. However, if the neu- 
trino chemical potential is sufficiently large [18|, the QCD 
phase transition does become first order. The required 
chemical potential does not violate nucleosynthesis con- 
straints, and if a sterile neutrino is the dark matter, we do 
actually expect a large neutrino chemical potential [19(. 

Pulsar timing experiments will reach unprecedented 
sensitivities in the next few years, and may open a new 
window on cosmology. The detection of a stochastic 
background with pulsar timing experiments could help 
to study the nature of the QCD phase transition, its du- 
ration, its strength an so on; by comparison with lattice 
calculations, this would allow us to determine the neu- 
trino chemical potential and other properties of the so 
elusive cosmological neutrino sector. Furthermore, the 
amplitude and peak frequency of the spectrum are sen- 
sitive to the expansion rate of the Universe at this tem- 
perature [20(| , which remains unconstrained to date. 

A first order QCD phase transition generating a GW 
background would also induce a stochastic background of 
magnetic fields, as studied in the past [561 ] . Furthermore, 
the GW background might have non-vanishing hclicity, 
which would be an interesting phenomenon to investigate 
by itself |3. 
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